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ABSTRACT 



Since cosmology is no longer "the data-starved science" , the problem of how to best 
analyze large data sets has recently received considerable attention, and Karhuncn- 
' Loeve eigenvalue methods have been applied to both galaxy redshift surveys and 

, Cosmic Microwave Background (CMB) maps. We present a comprehensive discussion 

(T) ■ of methods for estimating cosmological parameters from large data sets, which includes 

f^*) | the previously published techniques as special cases. We show that both the problem of 

^sO . estimating several parameters jointly and the problem of not knowing the parameters 

a priori can be readily solved by adding an extra singular value decomposition step. 
^ (— I It has recently been argued that the information content in a sky map from a next 

O generation CMB satellite is sufficient to measure key cosmological parameters (h, il, 

A, etc.) to an accuracy of a few percent or better — in principle. In practice, the 
data set is so large that both a brute force likelihood analysis and a direct expansion 
in signal-to-noise eigenmodes will be computationally unfeasible. We argue that it is 
likely that a Karhunen-Loeve approach can nonetheless measure the parameters with 
close to maximal accuracy, if preceded by an appropriate form of quadratic "pre- 
compression" . 

We also discuss practical issues regarding parameter estimation from present and 
future galaxy redshift surveys, and illustrate this with a generalized eigenmode analysis 
of the IRAS 1.2 Jy survey optimized for measuring (3 = tt°- 6 /b using redshift space 
distortions. 

Key words: Cosmology: theory, large scale structure of Universe, cosmic microwave 
background - Astronomical methods: data analysis, statistical 



1 INTRODUCTION 

The problem of analysis of large data sets is one that, un- 
til recently, has not been a major concern of cosmologists. 
Indeed, in some areas no data existed to be analyzed. In 
the last few years, this situation has rapidly changed. A 
highlight in this transition has been the discovery of fluc- 
tuations in the cosmic microwave background (CMB) by 
the Cosmic Background Explorer (COBE) satellite (Smoot 
et al. 1992). In its short lifetime, COBE produced such a 
large data set that a number of sophisticated data-analysis 
methods were developed specifically to tackle it. In addi- 
tion, the advent of large galaxy redshift surveys has created 
a field where the data sets increase by an order of magnitude 
in size in each generation. For instance, the surveys where 



the object selection was based on the Infrared Astronomy 
Satellite (IRAS) contain several thousand galaxies: ~ 2, 000 
(QDOT; Lawrence et al. 1996) ~ 5,000 (Berkeley 1.2Jy; 
Fisher et al. 1995) and ~ 15, 000 (PSC-z; Saunders et al. 
1996). The proposed next-generation surveys will have much 
larger numbers of objects — around 250,000 in the Anglo- 
Australian Telescope 2 degree Field galaxy redshift survey 
(Taylor 1995) and ~ 10 6 for the Sloan Digital Sky Survey 
(Gunn & Weinberg 1995). Similarly plate measuring ma- 
chines, such as the APM at Cambridge and SuperCOSMOS 
at the Royal Observatory, Edinburgh, can produce very large 
catalogues of objects, and numerical simulations of galaxy 
clustering are even now capable of producing so much data 
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that the analysis and storage of the information is in itself 
a challenge. 

A standard technique for estimating parameters from 
data is the brute force maximum likelihood method, which 
illustrates why people have been driven towards developing 
more sophisticated methods. For n data items {e.g., pixels 
in a CMB map, or Fourier amplitudes from a transformed 
galaxy distribution), the maximum likelihood method re- 
quires inversion of an n x n matrix for each set of parameter 
values considered — and this is for the simplest possible 
case where the probability distribution is Gaussian. Since a 
next-generation CMB satellite might produce a high resolu- 
tion sky map with ~ 10 7 pixels, and the CPU time required 
for an inversion scales as n 3 , a brute force likelihood analysis 
of this type of data set will hardly be feasible in the near 
future. 

Fortunately, it is often possible to greatly accelerate 
a likelihood analysis by first performing some appropriate 
form of data compression, by which the data set is substan- 
tially reduced in size while nonetheless retaining virtually all 
the relevant cosmological information. In this spirit, a large 
number of data-compression methods have been applied in 
the analysis of both CMB maps {e.g. Seljak & Bertschinger 
1993; Gorski 1994) and galaxy redshift surveys {e.g. Davis & 
Peebles 1983; Feldman et al. 1994; Heavens & Taylor 1995). 
A powerful prescription for how to do this optimally is the 
Karhunen-Loeve eigenvalue method (Karhunen 1947) , which 
has recently been applied to both CMB maps (Bond 1994; 
Bunn 1995; Bunn & Sugiyama 1995) and redshift surveys 
(Vogeley 1995; Vogeley & Szalay 1996). The goal of this pa- 
per is to review the more general framework in which these 
treatments belong, and to present some important general- 
izations that will facilitate the analysis of the next genera- 
tion of cosmological data sets. 

The rest of this paper is organized as follows. In Sec- 
tion [j| we review some useful information-theoretical results 
that tell us how well parameters can be estimated, and how 
to determine whether a given analysis method is good or 
bad. In Section^, we review the Karhunen-Loeve data com- 
pression method and present some useful generalizations. In 
Section |5| we illustrate the theory with various cosmology ap- 
plications, including the special case of the signal-to-noise 
eigenmode method. In Section Jf] we discuss limitations of 
method and possible ways of extending it to make the analy- 
sis feasible for huge data sets such as a 10 7 pixel future CMB 
map. Finally, our results are summarized and discussed in 
Section |(| 



2 HOW WELL CAN YOU DO WITHOUT 
DATA COMPRESSION? 

How accurately can we estimate model parameters from a 
given data set? This question was basically answered 60 
years ago (Fisher 1935), and we will now summarize the 
results, which are both simple and useful. 



2.1 The classical results 

Suppose for definiteness that our data set consists on n real 
numbers xi, X2, x n , which we arrange in an n-dimensional 



vector x. These numbers could for instance denote the mea- 
sured temperatures in the n pixels of a CMB sky map, the 
counts-in-cells of a galaxy redshift survey, the first n coef- 
ficients of a Fourier-Bessel expansion of an observed galaxy 
density field, or the number of gamma-ray bursts observed 
in n different flux bins. Before collecting the data, we think 
of a; as a random variable with some probability distribution 
L{x; 0), which depends in some known way on a vector of 
model parameters 



= {61,82, dn 



(1) 



Such model parameters might for instance be the spectral 
index of density fluctuations, the Hubble constant h, the cos- 
mic density parameter Q or the mean redshift of gamma-ray 
bursts. We will let ©0 denote the true parameter values and 
let refer to our estimate of 0. Since is some function 
of the data vector x, it too is a random variable. For it to be 
a good estimate, we would of course like it to be unbiased, 
i.e., 



{0} = ©o, 



(2) 



and give as small error bars as possible, i.e., minimize the 
standard deviations 



A9i = ((0?) - {9 t f) 



2\l/2 



(3) 



In statistics jargon, we want the BUE 8i, which stands for 
the "Best Unbiased Estimator". 

A key quantity in this context is the so-called Fisher 
information matrix, defined as 



d 2 c 

86,80, 



(4) 



where C = — \nL. Another key quantity is the maximum 
likelihood estimator, or ML-estimator for brevity, defined as 
the parameter vector @]y[L that maximizes the likelihood 
function L{x; 0). (Above we thought of L{x; 0) as a prob- 
ability distribution, a function of x. However, when limiting 
ourselves to a given data set x and thinking of L{x; 0) as 
a function of the parameters 0, it is customarily called the 
likelihood function.) 

Using this notation, a number of powerful theorems 
have been proven (see e.g. Kenney & Keeping 1951 and 
Kendall & Stuart 1969 for details): 

(i) For any unbiased estimator, Adi > l/\/Fa. 

(ii) If there is a BUE 0, then it is the ML-estimator or a 
function thereof. 

(iii) The ML-estimator is asymptotically BUE. 

The first of these theorems, known as the Cramer-Rao in- 
equality, thus places a firm lower limit on the error bars 
that one can attain, regardless of which method one is us- 
ing to estimate the parameters from the data. This is the 
minimum error bar attainable on 0i if all the other parame- 
ters are known. If the other parameters are estimated from 
the data as well, the minimum standard deviation rises to 
AS, > (F-^ 1 / 2 . 

The second theorem shows that maximum-likelihood 
(ML) estimates have quite a special status: if there is a best 
method, then the ML-method is the one. Finally, the third 
result basically tells us that in the limit of a very large data 
set, the ML-estimate for all practical purposes is the best 
estimate, the one that for which the Cramer-Rao inequality 
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becomes an equality. It is these nice properties that have 
made ML-estimators so popular. 



2.2 The Fisher information matrix 

Although the proof of the Cramer-Rao inequality is rather 
lengthy, it is quite easy to acquire some intuition for why the 
Fisher information matrix has the form that it does. This is 
the purpose of the present section. 

Let us Taylor expand C around the ML-estimate 0. 
By definition, all the first derivatives dC/d6i will vanish at 
the ML-point, since the likelihood function has its maxi- 
mum there, so the local behavior will be dominated by the 
quadratic terms. Since L = exp[— £,], we thus see that the 
likelihood function will be approximately Gaussian near the 
ML-point. If the error bars are quite small, L usually drops 
sharply before third order terms have become important, so 
that this Gaussian is a good approximation to L everywhere. 
Interpreting L as a Bayesian probability distribution in pa- 
rameter space, the covariance matrix T is thus given simply 
by the second derivatives at the ML-point, as the inverse of 
the Hessian matrix: 



(T- 



d 2 c 



(•») 



Note that the Fisher information matrix F is simply the 
expectation value of this quantity at the point = 0o 
(which coincides with the ML-point on average if the ML- 
estimate is unbiased). It is basically a measure of how fast 
(on average) the likelihood function falls off around the ML- 
point, i.e., a measure of the width and shape of the peak. 
From this discussion, it is also clear that we can use its 
inverse, i* 1-1 , as an estimate of the covariance matrix 



T = {00 t )-{0){0) t 



(6) 



(ii) 



Since C is a symmetric matrix for all values of the parame- 
ters, it is easy to see that all the derivatives C,i, C etc., 
will also be symmetric matrices. Using the matrix identities 
(C- 1 ),i = -C~^C,iC- x and (lnC%= C^C,,, we thus 
obtain 

2C, r = Tr [C^C,* -C~ x C,i C X D + C X B,, ] . (12) 

When evaluating C and fi at the true parameter values, we 
have (x) = fj, and (xx l ) = C + fAfl 1 , which gives 

r (d) = c, 

=0, (13) 
{D,ij) =/x,i/x,5+/x,j/x,'. 

Using this and equation (|l2j), we obtain {C,i ) = 0. In other 
words, the ML-estimate is correct on average in the sense 
that the average slope of the likelihood function is zero at the 
point corresponding to the true parameter values. Applying 
the chain rule to equation (|l^), we obtain 

1£,,ij — Tr[ C C ,i C C ,j ~\~C C ,ij 

+ C- 1 (C, i C- 1 C, :i +C, j C- 1 C,i)C- 1 D 

- C- 1 (C, i C- 1 D, :i + C, j C- 1 D,i) 

- C- 1 C, lJ C- 1 D + C- 1 D, lJ ]. (14) 

Substituting this and equation ([l3|) into equation (Q) and 
using the trace identity Tr[AB] = Tr[B A] , many terms drop 
out and the Fisher information matrix reduces to simply 



(15) 



of our parameter estimates when we use the ML-method. 



where we have defined the matrices A* = C~ 1 C,i — (In C),j 
and Mij = (D,ij ) = M» fl,j +H,j fJ>,t- This old and well- 
known result is also derived by Bunn (1995) and Vogeley & 
Szalay (1996). 



2.3 The Gaussian case 

Let us now explicitly compute the Fisher information ma- 
trix for the case when the probability distribution is Gaus- 
sian, i.e., where (dropping an irrelevant additive constant 
nln[27r]) 



2C = IndetC + (x - fi)C' 1 {x - fx)* 



where in general both the mean vector fi and the covariance 
matrix 



c = <(* - n)( x - n)*) 



(8) 



depend on the model parameters 0. Although vastly sim- 
pler than the most general situation, the Gaussian case is 
nonetheless general enough to be applicable to a wide variety 
of problems in cosmology. Defining the data matrix 



D = (x - n)(x - 



(9) 



and using the well-known matrix identity In det C = Tr In C, 
we can re-write equation (Q) as 



2C = Tr [in C + C^D] 



(10) 



We will use the standard comma notation for derivatives, 
where for instance 



2.4 An example: parameter estimation with a 
future CMB experiment 

Let us illustrate the above results with a simple example 
where fi = and C is diagonal. Suppose a next generation 
CMB satellite generates a high-resolution all-sky map of the 
CMB fluctuations. Let us for the moment ignore the com- 
plication of the galactic plane, and expand the temperature 
distribution in spherical harmonics ai m . We define our data 
vector x as the set of these coefficients from £ = 2 up to some 
cutoff £ m ax, i.e., n = (£max + l) 2 -4 and Xi = ai, where we 
have combined I and m into a single index i = 1,2, 3... as 
i — £ 2 + £ + m + 1. With the standard assumption that the 
CMB fluctuations are isotropic, we have 

( fi =0, 



C ij — &ij 



(16) 



Here Ci denotes the angular power spectrum of the CMB, 
and the second term incorporates the effects of instrumen- 
tal noise and beam smearing (Knox 1995, Tegmark & Efs- 
tathiou 1996). iV is the number of pixels in the sky map, a 
is the r.m.s. pixel noise, 



= FWHM/V81n2 w 0.425 FWHM, 



(17) 
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Figure 1. The derivatives of the CDM power spectrum with 
respect to various parameters. 



and FWHM is the full-width-half-max beam width. Assum- 
ing that the CMB and noise fluctuations are Gaussian, equa- 
tion ([y| gives the Fisher information matrix 



'max 



47TO- 2 9^1(1+1) 

N 6 



( dC t 



dCt 



This handy expression, also derived by Jungman et al. 
(1996a), tells us that the crucial functions which determine 
the attainable accuracy are the derivatives of the power spec- 
trum with respect to the various parameters. Examples of 
such derivatives are shown in Figure |l|. For instance, the 
derivative with respect to the reionization optical depth, 
dCi/dr, is shaped as — Ct for £ 3> 10, since earlier reion- 
ization suppresses all these multipoles by the same factor 
e -2r rp^g derivative with respect to the quadrupole nor- 
malization, dCe/dQ, of course has the same shape as the 
power spectrum itself. 

Equation ( |l8| ) has a simple geometric interpretation. We 
can think of the m functions dCe/dOi as vectors in a Hilbert 
space of dimension (imax— 1), and think of the Fisher matrix 
component Fij as simply the dot product of the vectors 
dCi/dOi and dCi/d9i. This dot product gives a rather small 
weight to low ^-values, essentially because there are fewer m- 
modes there and correspondingly larger cosmic variance. In 
addition, the weights are seen to be exponentially suppressed 
for large I, where beam smearing causes the effect of pixel 
noise to explode. If the parameter dependence of Ce was 
such that all n vectors dCi/d6i were orthogonal under this 
dot product, then F and the parameter covariance matrix 
T — F^ 1 would be diagonal, and the errors in the estimates 
of the different parameters would be uncorrelated. The more 
similarly shaped two curves in Figure |l| are, the harder it 
will be to break the degeneracy between the corresponding 
two parameters. In the extreme case where two curves have 



identical shape (are proportional to each other), the Fisher 
matrix becomes singular and the resulting error bars on the 
two parameters become infinite. It is therefore interesting to 
diagonalize the matrix T (or equivalently, its inverse F). The 
eigenvectors will tell us which m parameter combinations 
can be independently estimated from the CMB data, and 
the corresponding eigenvalues tell us how accurately this 
can be done. 

It has been pointed out (Bond et al. 1994) that there 
will be a considerable parameter degeneracy if data is only 
available up to around the first Doppler peak. This is clearly 
illustrated in Figure most of the curves lack strong fea- 
tures in this range, so some of them can be well approxi- 
mated by linear combinations of the others. If a CMB ex- 
periment has high enough resolution to measure the power 
out to I ~ 10 3 , however, this degeneracy is broken. Jungman 
et al. (1996b) compute the Fisher matrix for a model with 
the eleven unknown parameters 



= (fi, U b h 2 , h, A, n s , a, n T , T/S, r, Q,N V ); 



(19) 



the density parameter, the baryon density, the Hubble pa- 
rameter, the cosmological constant, the spectral index of 
scalar fluctuations, the "running" of this index (a measure 
of the deviation from power law behavior), the spectral in- 
dex of tensor fluctuations, the quadrupole tensor-to-scalar 
ratio, the optical depth due to reionization, and the number 
of light neutrino species, respectively. They find that even 
when estimating all these parameters simultaneously, the re- 
sulting error bars on some of the key parameters are of the 
order of a few percent or better for experiments with a res- 
olution good enough to probe the first few Doppler peaks, 
n&ren the abundance of hot dark matter can be constrained 
in this fashion (Dodelson et al. 1996). An up-to-date discus- 
sion of which additional parameters can be constrained when 
using large-scale-structure data as well is given by Liddle et 
al. (1996). 

In should be stressed that the formalism above only tells 
us what the attainable accuracy is if the truth is somewhere 
near the point in parameter space at which we compute the 
power derivatives. Figure |l| corresponds to standard COBE- 
normalized CDM, i.e., to 



0o = (1, 0.015, 0.5, 0, 1, 0, 0, 0, 0, 20^K, 3). 



(20) 



The worst scenario imaginable would probably be extremely 
early reionization, since r ^> 1 would eliminate almost all 
small-scale power and introduce a severe degeneracy by mak- 
ing all the power derivatives almost indistinguishable for 
£ ^> 10, the region where they differ the most in Figure |l| 

Needless to say, accurate parameter estimation also re- 
quires that we can compute the theoretical power spectrum 
accurately. It has been argued (Hu et al. 1995) that this can 
be done to better than 1%, by accurately modeling various 
weak physical effects such as Helium recombination and po- 
larization. We will return to other real-world issues, such 
as foreground contamination and incomplete sky coverage, 
further on. 
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3 OPTIMAL DATA COMPRESSION: 
KARHUNEN-LOEVE METHODS 

Above we saw how small error bars one could obtain when 
using all the information in the data set. We also saw that 
for large data sets, these minimal error bars could be ap- 
proximately attained by making a likelihood analysis using 
the entire data set. However, there are of course many situ- 
ations where this approach is not feasible even in the simple 
Gaussian case. 



3.1 The need for data compression 

If we wish to estimate m parameters from n data points, 
this would involve inverting annxn matrix at a dense grid 
of points in the m-dimensional parameter space. The CPU 
time required to invert such a non-sparse matrix scales as 
n 3 , and the number of grid points grows exponentially with 
m, so there are obviously limits as to what can be done in 
practice. 

For instance, the "brute force" COBE analysis of 
Tegmark & Bunn (1995) had n = 4038 and m = 2, and 
involved inverting 4038 x 4038 matrices at a couple of hun- 
dred grid points. A future high-resolution all-sky CMB map 
might contain of order n — 10 7 pixels, and direct numerical 
inversion of matrices this large is clearly unfeasible at the 
present time. Also, to numerically map out the likelihood 
function in the 11-dimensional parameter space explored by 
Jungman et al. (1996b) with reasonable resolution would 
entail evaluating it at such a large number of grid points, 
that it would be desirable to make the inversion at each 
grid point as fast as possible. 

Likewise, the spherical harmonic analysis of the 1.2Jy 
redshift survey (Heavens & Taylor, 1995, Ballinger Heav- 
ens & Taylor 1996), and the PSC-z survey (Tadros et al., 
in preparation), required a likelihood analysis of some 1208 
modes to find the redshift distortion parameter and a non- 
parametric stepwise estimate of the power spectrum. 

Because of this, it is common to carry out some form of 
data compression whereby the n data points are reduced to 
some smaller set of n' numbers before the likelihood analysis 
is done. Since the time needed for a matrix inversion scales as 
n 3 , even a modest compression factor such as n/n' = 5 can 
accelerate the parameter estimation by more than a factor 
of 100. 



3.2 The optimization problem 

In this section, we will discuss the special case of linear data 
compression, which has as an important special case the so- 
called signal-to-noise eigenmode method. We will return to 
non-linear data compression methods below, in Section 

The most general linear data compression can clearly 
be written as 



y = Bx, 



(21) 



where the n'-dimensional vector y is the new compressed 
data set and B is some arbitrary n' x n matrix. Thus the 
new data set has 



(y) 

(yy t )-{y)(y) t 



= B\x, 
= BCB\ 



(22) 



and substituting this into equation ( |15[ ) we find that the new 
Fisher information matrix F is given by 



= 2 Tr ' 



(BCB t )~ 1 (BC,i B t )(BCB t )~ 1 (BC,j B*) 



+ (BCB )~ (BMijB )]. (23) 

If n = n' and B is invertible, then the B-matrices cancel in 
this expression, leaving us with simply 

1 



^ = -Tr[B" 



\AiAj + C- 1 M li )B t ] = Fy 



(24) 



In other words, B just makes a similarity transform of the 
matrix within the trace in equation (|l5|), leaving the Fisher 
information matrix F unchanged. This reflects the fact that 
when B is invertible, no information is lost or gained, so 
that the error bars on a measured parameter are unchanged. 
For instance, replacing a galaxy-cut COBE map consisting 
of say 3600 pixels by its spherical harmonic coefficients up 
to £ = 59 (there are 3600 such coefficients) would be an 
invertible transformation, thus making no difference what- 
soever as to the attainable error bars. Likewise, expansion 
in galaxy-cut spherical harmonics gives exactly the same re- 
sult as expansion in an orthonormal basis for the galaxy-cut 
spherical harmonics (as in Gorski 1994), since the mapping 
from the non-orthonormal basis to the orthonormal basis is 
invertible. This result is obvious if we think in terms of in- 
formation: if the mapping from x to y is invertible, then y 
must clearly contain exactly the same information that x 
does, since we can reconstruct x by knowing y. 

Rather, the interesting question is what happens when 
n < n and B is not invertible. Each row vector of B spec- 
ifies one of the numbers in the new data set (as some lin- 
ear combination the original data x), and we wish to know 
which such row vectors capture the most information about 
the parameters 0. If B has merely a single row, then we can 
write B = b* for some vector b, and the diagonal (i — j) 
part of equation (0) reduces to 



- _ 1 [ b'Cib 
ll ~ 2 { b'Cb 



+ 



(b'Cb) 



(25) 



Let us now focus on the problem of estimating a single pa- 
rameter 6i when all other parameters are known — we will 



return to the multi-parameter case below, in Section 3.6. 

1/2 

Since the error bar A0i w F u in this case, we want to find 
the b that maximizes the right-hand side of equation (pHf), 
Although this is a non-linear problem in general, there are 
two simple special cases which between them cover many of 
the situations that appear in cosmology applications: 

• The case where the mean is known: fi,i = 

• The case where the covariance matrix is known: C,i = 

Below we show first how these two cases can be solved sepa- 
rately, then how the most general case can be for all practical 
purposes solved by combining these two solutions. 



3.3 When the mean is known 

When u, is independent of 6i, the second term in equa- 
tion (B3) vanishes, and we simply wish to maximize the 
quantity 

1/2 _ \b f C,j b\ 



(2Fii) 



b'Cb 



(26) 
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Since the mean /J, does not depend on the parameters (i.e., 
it is known), let us for simplicity redefine our data by 
x I— > x — fi, so that (a;) = 0. Note that whereas the de- 
nominator in equation ( ptj ) is positive (since C is a co- 
variance matrix and thus positive definite), the expression 
b l C ,i b in the numerator might be negative, since C,j is 
not necessarily positive definite. We therefore want to make 
(6'C,j &)/(&* C7f>) either as large as possible or as small as 
possible, depending on its sign. Regardless of the sign, we 
thus seek the b for which this ratio takes an extremum, i.e., 
we want the derivatives with respect to the components of b 
to vanish. Since this ratio is clearly unchanged if we multiply 
b by a constant, we can without loss of generality choose to 
normalize b so that the denominator equals unity. We thus 
seek an extremum of b l C,i b subject to the constraint that 



VCb 



(27) 



This optimization problem is readily solved by introducing 
a Lagrange multiplier A and extremizing the Lagrangean 

b'C^b- Xb'Cb. (28) 

Varying b and setting the result equal to zero, we obtain the 
generalized eigenvalue problem 



C,ib= XCb. 



(29) 



Since C is symmetric and positive definite, we can Cholesky 
decompose it (e.g., Press et al. 1992) as C = LL l for some 
invertible matrix L. Multiplying equation ( pg| ) by L^ 1 from 
the left, we can rewrite it as 



(L-'C, 



L'^lSb) = X(L%), 



(30) 



which we recognize as an ordinary eigenvalue problem for 
the symmetric matrix (L~ 1 C,i L~*), whose solution will be 
n orthogonal eigenvectors (L l b k ) and corresponding eigen- 
values Xk, k — 1,2, n. Let us sort them so that 



|Ai| > |A a | > ... > |A„|. 



(31) 



Let us choose the k th row of B to be the row vector b\ , so 
that the compressed data set is given by y k = b\x. The or- 
thogonality property (L bk) ■ (L • bk') oc 5 kk i in combination 
with our chosen normalization of equation (^) then tells us 
that our compressed data satisfies 



(VkVk>) = {(b t k x)(x t b k ,)} = b\Cb k , = b^LVb^, 



,(32) 



i.e., (yy 1 ) = /. The compressed data values y thus have the 
nice property that they are what is known as statistically or- 
thonormal, i.e., they are all uncorrelated and all have unit 
variance. Since we are assuming that everything is Gaussian, 
this also implies that they are statistically independent — in 
fact, y is merely a vector of independent unit Gaussian ran- 
dom variables. In other words, knowledge of y k gives us no 
information at all about the other y- values. This means that 
the entire information content of the initial data set x has 
been portioned out in disjoint pieces in the new compressed 
data set, where yi contains the most information about 9i, 
y2 is the runner up (once the information content of yi has 
been removed, j/2 contains the most information about 9i), 
and so on. 

If we use all the n vectors b k as rows in B, then B will be 
invertible, and we clearly have not achieved any compression 
at all. However, once this matrix B has been found, the 
compression prescription is obvious: if we want a compressed 



data set of size n' < n, then we simply throw away all but 
the first n' rows of B. It is straightforward to show (see 
e.g. Therrien 1992; Vogeley & Szalay 1996) that if we fix an 
integer n' and then attempt to minimize A9i, we will arrive 
at exactly the same B as with our prescription above (or 
this B multiplied by some invertible matrix, which as we 
saw will leave the information content unchanged). 

How does the error bar A8i depend on the choice of n'? 
Equation (|29|) implies that 

(33) 



C^B* = CB'A, 



where Aij = SijXi, i.e., a diagonal matrix containing all the 
eigenvalues. Since equation (B2I) implies that 

BCB 1 = I, 



it follows from equation (|33|) that 
BC.i B f = A. 



(34) 



(35) 



Therefore the diagonal part of equation (|23|) reduces to sim- 
ply 



2F U = Tt[{(BC B*)^ 1 (BC ,i B')} 2 ] = TrA 2 



2 ^. (36) 



It is thus convenient to plot X\ as a function of the mode 
number k, as we have done in Figure |^ for the specific case of 
measuring the redshift distortion parameter from the 1.2Jy 
redshift survey, to see at what k the information content 
starts petering out. Alternatively, one can plot A9i as a func- 
tion of n' as in Figure [|and see when the curve flattens out, 
i.e., when the computational cost of including more modes 
begins to yield diminishing returns. If one feels uncomfort- 
able about choosing n' by "eyeballing" , one can obtain a 
more quantitative criterion by noting that using only a sin- 
gle mode b k would give A9i = l/|Afc|. Thus the quantity 
0i|Afc| can be thought of as the signal-to- noise ratio for the 
fc th mode, and one can choose to keep only those modes 
where this ratio exceed unity (or some other constant such 
as 0.2). 

3.4 When the covariance matrix is known 

When C is independent of 9i , the first term in equation ( |25[ ) 
vanishes, and we simply wish to maximize the quantity 



p.. — 



b l M lt b 



(37) 



VCb 

Introducing a Lagrange multiplier A and proceeding exactly 
as in the previous section, this is equivalent to the eigenvalue 
problem 



(L _1 Af = X(L*b). 



(38) 



However, since the matrix Ma = 2/x,; fi,j is merely of rank 
one, this problem is much simpler than that in the previous 
section. Since the left hand side is 



2(L- 1 f J,, i )(fX,tb) oc (L 1 n„ 



(39) 



it points in the direction (L fl,i ) regardless what the vec- 
tor b is, so the only non-trivial eigenvector of (L~ l MaL - *) 
is 



(L t b 1 ) = (L- 1 n, i ), 



(40) 
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with a corresponding eigenvalue 
Ai = 2\L~ 1 fi, i 2 = Tr[CMu]. 



(41) 



All other eigenvalues vanish. In other words, b\ = C~ fj,,i, 
and the new compressed data set consists of just one single 
number, 



(42) 



which contains just as much information about the param- 
eter 9i as the entire data set x did. Another way to see this 
is to compute the Fisher information before and after the 
compression, and note that 



Fa — Fa — jLt,* C 1 fi, 



(43) 



In other words, all the information about 6i in x is has 
been distilled into the number y\. Thus the ML-estimate of 
8i is just some function of y\. For the special case where 
fi depends linearly on the parameter, i.e., fi = fi,i 9i where 
is a known constant vector, this function becomes trivial: 
the expectation value of equation (^) reduces to 

<y 1 ) = (At,*C- 1 At, ! )e i , (44) 

which means that on average, up to a known constant 
(/l,jC~ the compressed data is just the desired pa- 

rameter itself. 



3.5 When neither is known — the general case 

In the most general case, both fi and C depend on the 
parameters 0. This happens for instance in the CMB case 
when the components of x are chosen to be estimates of the 
power spectrum Ct, as in Section |B| below. 

Clearly, the more ways in which the probability distri- 
bution for x depends on the parameters, the more accurately 
we can estimate them. In other words, when both fM and C 
depend on 0, we expect to do even better than in the two 
cases discussed above. Finding the vectors b that extremize 
equation (^H|) is quite a difficult numerical problem. For- 
tunately, this is complet ely unnecessary. Solving the eigen- 



value problem in Section 3.3 gives us n compression vectors 
capturing the bulk of the cosmological information coming 
from the first term in equation (t25J). Section i.4 gives us 
one additional vector C _1 that captures all the cosmo- 
logical information coming from the second term. In other 
words, if we simply use all these n" + 1 compression vectors 
together, we have for all practical purposes solved our prob- 
lem. Since a direct numerical attack on equation ( p5|) could 
never reduce these n" + 1 vectors to fewer than n" without 
widening the resulting error bars, the time savings in the 
ensuing likelihood analysis would at most be a completely 
negligible factor [(n" + l)/n"] 3 ~ 1 + 3/n" compared the 
simple prescription we have given here. 



3.6 Estimating several parameters at once 

The Karhunen-Loeve method described above was designed 
to condense the bulk of the information about a single 
parameter into as few numbers yt as possible. Although 
this particular problem does occasionally arise, most cos- 
mology applications are more complicated, involving mod- 
els that contain more than one unknown parameter. The 



previous applications of Karhunen-Loeve methods to CMB 
maps (Bond 1994; Bunn & Sugiyama 1995; Bunn 1995) and 
galaxy surveys (Vogeley 1995; Vogeley & Szalay 1996), have 
involved taking a rather minimalistic approach to this issue, 
by optimizing the eigenmodes for one particular parameter 
(typically the overall normalization of the fluctuations) and 
assuming that these eigenmodes will capture most of the 
relevant information about the other parameters as well. 

If we want to estimate several parameters from a data 
set simultaneously, then how should we best compress the 
data? As we saw in Section |^, the covariance matrix T of 
our parameter estimates is approximately -F" 1 , so for the 
multivariate case, we have 



A6i = (F 



-l\l/2 



(45) 



instead of A8i — F it . One approach to the problem would 
thus be trying to minimize the diagonal elements of F~ 
rather than, as above, trying to minimize the diagonal ele- 
ments of F. This is, however, quite a cumbersome optimiza- 
tion problem. 

Fortunately, there is an alternative solution to the prob- 
lem which is easy to implement in practice and in addition is 
easy to understand intuitively, in terms of information con- 
tent. Suppose that we repeat the standard KL procedure 
described above m times, once for each parameter 8i, and 
let n'i denote the number of modes that we choose to use 
when estimating the i th parameter. We then end up with 
m different compressed data sets, such that the n'i numbers 
in the i th set contain essentially all the information that 
there was about 9i. This means that the union y of all these 
compressed data sets, consisting of 



J2< 



(46) 



numbers, retains basically all the information from x that is 
relevant for our joint estimation of all the parameters. The 
problem is of course that there might be plenty of redun- 
dancy in y. Indeed, if we typically compressed by say a fac- 
tor n/n" ~ 10 and had 11 parameters as in the CMB model 
of Jungman et al. (1996b), then we would have achieved a 
counterproductive "anti-compression" with n" > n. 

How can we remove the unnecessary "pork" from y? 
The n'i eigenvectors selected via the KL compression method 
span a certain subspace of the n-dimensional space in which 
the data vector x lives, and we can think of the data com- 
pression step as simply projecting the vector x onto this 
subspace. For each parameter 8i , we found such a subspace, 
and the redundancy stems from the fact that many of these 
subspaces overlap with one another, or point in very simi- 
lar directions. In order to compress efficiently and yet retain 
almost all the the information about all the parameters, we 
want to find a small set of vectors that span as much as pos- 
sible of the union of all the subspaces. As described in detail 
in e.g. Press et al. (1992), this is readily accomplished by 
making a Singular Value Decomposition (SVD) and throw- 
ing away all vectors corresponding to very small singular 
values. Let us define Bi = BA, where B is the KL com- 
pression matrix optimized for estimating the i parameter. 
Here we multiply the row vectors by their corresponding 
eigenvalues so that better vectors will receive more weight 
in what follows. Combining all the row vectors of the trans- 
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posed compression matrices B\ into a single matrix, SVD 
allows us to factor this matrix as 

(B\ ... B* m ) =UAV t , (47) 

where U t U — I, V V = I and the matrix A is diagonal. 
One customarily writes A = diag{Ai}, and refers to the di- 
agonal elements A* as the singular values. These are basically 
a generalization of eigenvalues to non-square matrices. We 
define our final compression matrix as 

B = U l . (48) 

The columns of U with corresponding non-zero singular val- 
ues Xi form an orthonormal set of basis vectors that span the 
same space as all the initial compression vectors together. 
The column vectors of U with vanishing singular values form 
an orthonormal basis for the null-space of UAV 1 , repre- 
senting redundant information. Similarly, the vectors corre- 
sponding to singular values near zero capture information 
that is almost entirely redundant. By making a plot similar 
to Figure |^, showing \k as a function of k, one decides on 
the final number n' < n" of row vectors in B to keep. Once 
n' has been fixed, one can of course (if one prefers^) make 
the n' numbers in the new compressed data set statistically 
orthogonal as before by Cholesky decomposing their covari- 
ance matrix BCB 1 — LL l and replacing z = Bx by L~ f z. 
In summary, we have found that when we wish to estimate 
more than one parameter from the data, we can obtain a 
close to optimal data compression with the following steps: 

(i) Compute the KL eigenmodes that contain the bulk of 
the information about the first parameter. 

(ii) Repeat this procedure separately for all other param- 
eters. 

(iii) Arrange all the resulting vectors, multiplied by their 
eigenvalues, as the rows of B. 

(iv) Make an SVD of B and throw away all rows corre- 
sponding to very small singular values. 

3.7 The problem of not knowing the parameters a 
priori 

The KL-approach has sometimes been criticized for not be- 
ing model-independent, in the sense that one needs to make 
an a priori guess as to what the true parameter values 0q 
are in order to be able to compute the compression matrix 
B. Although there is no evidence that a bad guess as to 0q 
biases the final results (see e.g. Bunn 1995 for a detailed 
discussion of these issues, including numerical examples), a 
poor guess of 0q will of course lead to a data compression 
that is no longer optimal. In other words, one would expect 
to still get an unbiased answer, but perhaps with larger error 
bars than necessary. 

In practice, this loss of efficiency is likely to be marginal. 
For instance, as a test, Bunn (1995) performed a KL-analysis 
of a COBE CMB map with n ~ 1 assuming a blatantly 
erroneous spectral index n — 2 to compute B (n — 2 is 
ruled out at about 3<r a posteriori), and compared this with 

* The only merit of the statistical orthogonality is that it will 
make (yy 1 ) more close to diagonal near the fiducial parameter 
values 0, which might conceivably speed up the matrix inversion 
if an iterative method is used. 



the results obtained when assuming n = 1. The error bars 
where found to change only marginally. 

If one nonetheless wishes to do something about this 
efficiency problem, it is fortunately quite easy in practice. 
The simplest approach is an iterative scheme, whereby the 
KL-analysis is carried out twice, using the ML-estimate of 
from the first run as the fiducial "true" value in the second 
run. A more rigorous approach is to compute compression 
vectors b for a number of different assumptions about 0q 
that include the extreme corners of parameter space, and 
then combine all these vectors into a single compression ma- 
trix B by singular- value decomposition exactly as described 
in the previous section. 



4 COSMOLOGY APPLICATIONS 

In this section, we illustrate the theoretical discussion above 
with a number of cosmology examples, and show how the 
recently published work on CMB maps and galaxy surveys 
fits into the general KL-scheme as special cases. We will see 
that the typical compression factor is about 10 in both an 
IRAS example and a COBE example. 

4.1 A large-scale-structure example: 
redshift-space distortions 

Here we apply a KL-analysis to the problem of redshift-space 
distortions in the 1.2Jy IRAS galaxy survey (Fisher et al. 
1995). This problem has previously been analyzed in detail 
by Heavens & Taylor (1995, hereafter HT95) who used a 
spherical harmonic analysis. Here we repeat their analysis, 
but include a KL data-compression step to investigate how 
many modes are needed for an accurate determination of 
the redshift distortion parameter 



where b denotes the conventional linear bias factor, the ratio 
between the fluctuations in luminous matter and the total 
matter fluctuations. As was first shown by Kaiser (1987), the 
peculiar motions of galaxies induces a characteristic radial 
distortion in the apparent clustering pattern that depends 
only on this parameter /3. 

As the initial data vector x, we use the coefficients 
obtained by expanding the observed galaxy distribution in 
combinations of spherical harmonics and spherical Bessel 
functions as described in HT95, with the known mean val- 
ues subtracted off so that fi — (x) = 0. HT95 show that the 
covariance matrix is given by 

C* M „ = A™ + i ^($ n(l + pv aii ){& av +pV av )P{k a ),{50) 

a 

where the indices fi, v and a run over the above-mentioned 
modes, and P(k) is the power spectrum. The matrices <P 
and V encapsulate the effects of finite survey volume con- 
volution and redshift distortions, respectively. The matrix 
7l sn contains the shot-noise contribution. All three matrices 
are independent of /3. For our illustrative example, we as- 
sume a parametrised form of the power spectrum suggested 
by Peacock & Dodds (1994), which means that /3 is the only 
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Mode Number 
Figure 2. KL-eigenvalues A = l/A/3. 

unknown parameter. Hence m = 1 and & = 9i = f3. Differ- 
entiating equation (pp]), we obtain 

c,i = H = i(W + #py') + (vpv')/3, (5i) 

where P a /3 = P(k a )S a p (i.e., P = diag{P(fc a }). We assume 
an a priori value /3 = 1 to evaluate this matrix. Using n = 
1208 modes in the initial data set x as in HT95, we obtain 
the eigenvalues Xk shown in Figure []. The resulting error 
bar 



A/3 



E*2 



-1/2 



(52) 



is plotted as a function of n' (the number of modes used) in 
Figure and is seen to level out at around n' = 100. This 
means that although HT95 used ~ 10 3 modes to estimate f3, 
almost all the information is actually contained in the best 
~ 10 2 linear combinations of these modes. In other words, 
we can obtain basically identical results to those found in 
HT95 by using a compressed data set only 10% of the orig- 
inal size. Since the matrix inversion time scales as n' 3 , this 
compression factor of 10 thus allows the inversion to be car- 
ried out 1000 times faster. 



4.2 A special case: the signal-to-noise eigenmode 
method 

The above-mentioned example was rather generic in the 
sense that C depended on in a nonlinear way (in this 
case, C depended quadratically on /3). However, the special 
case of the KL-method where the parameter dependence is 
linear (or rather afEne) is so common and important that it 
has acquired a special name: the signal-to-noise eigenmode 
method. This is the special case where fi = 0, we have only 
one unknown parameter 9, and the covariance matrix can 
be written in the form 



S8 + N, 



(53) 



where the known matrices S and N are independent of 9. 
For reasons that will become clear from the examples below, 
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Figure 3. Error bar on beta (3 as a function of n' , the number 
of eigenmodes used. 



S and N are normally referred to as the signal and noise 
matrices, respectively, and they are normally both positive 
definite. Since dC/d9 — S, equation ( pj| ) gives the general- 
ized eigenvalue problem 



Sb = \(S + N)b. 



(54) 



By Cholesky decomposing the noise matrix as N = LI/. 
this is readily rearranged into the ordinary eigenvalue prob- 
lem 



(L~ X S , L~ t ){L t b) = \'{L*b), 



(55) 



where A' = A/(l — A). Since the matrix to be diagonal- 
ized, (L'^SL-*), is loosely speaking (N- 1/2 SN- 1/2 ), a 
type of signal-to-noise ratio, its eigenvectors b are usually 
referred to as the signal-to-noise eigenmodes. Historically, 
this terminology probably arose because one was analyzing 
one-dimensional time series with white noise, where AT oc I, 
so that (L^ 1 SL~ f ) — SN -1 was the signal-to-noise ratio 
even in the strict sense. Vogeley & Szalay (1996) refer to 
the change of variables x t— > L~ x x as prewhitening, since 
it transforms S + N into (L~ l SL~ l ) + I, i.e., transforms 
the noise matrix into the white noise matrix, the identity 
matrix. 

Equation (^2|) showed that in the general KL-case, the 
compressed data set was statistically orthonormal, {yy 1 } = 
I. In the S/N-case, the compressed data has an additional 
useful property: both the signal and noise contributions to 
y are uncorrelated separately. It is easy to show that 



bkSby = A' 'b\Nby oc 8^ 



(56) 



which means that the covariance matrix (yy ) will remain 
diagonal even if we have misestimated the true signal-to- 
noise ratio. In other words, if the sole purpose of our likeli- 
hood analysis is to estimate the overall normalization 9 from 
equation (|H^), we only need to invert diagonal matrices. 

The signal-to-noise method arises in a wide variety of 
contexts. For instance, it is a special case not only of the KL- 
method, but also of the power spectrum estimation method 
of Tegmark (1996) , corresponding to the case were the width 
of the window functions is ignored. 



10 M. Tegmark, A. N. Taylor, & A. F. Heavens 



4-2.1 Signal-to-noise analysis of CMB maps 

The signal-to-noise eigenmode method was introduced into 
CMB analysis by Bond (1994), who applied it to sky maps 
from both the COBE and FIRS experiments. It has also 
been applied to the COBE data by Bunn|], and used it to 
constrain a wide variety of cosmological models (Bunn 1995; 
Bunn & Sugiyama 1995; Bunn, Scott & White 1995; Hu, 
Bunn & Sugiyama 1995; White & Bunn 1995; Yamamoto & 
Bunn 1996). Here the uncompressed data set consists of the 
CMB temperatures in n pixels from a sky map (for instance, 
n = 4038 or 4016 for a COBE map with a 20° galactic cut, 
depending on the pixelization scheme used). If foreground 
contamination is negligible (see e.g. Tegmark & Efstathiou 
1996 for a recent review) and the pixel noise is uncorrelated 
(as it is for COBE to an excellent approximation — see 
Lineweaver et al. 1994) one has 



H =0, 
Nij = Sijo-j, 



(57) 



where At is a unit vector pointing the direction of the i th 
pixel, ai is the r.m.s. noise in this pixel, Pe are the Legen- 
dre polynomials and We is the experimental beam function. 
This expression for S should only be used if one marginal- 
izes over the monopole and dipole in the ensuing likelihood 
analysis — otherwise S should be corrected as described in 
Tegmark & Bunn (1995). In all the above-mentioned appli- 
cations, the compression was optimized with respect to the 
overall normalization of the power spectrum (say 9 = C2), 
so the matrix S was independent of 9. It has been found 
(Bunn & Sugiyama 1995) that a compression factor of about 
10, down to n' ~ 400 modes, can be attained without signifi- 
cant loss of information. This is about the same compression 
factor that we obtained in our redshift distortion example 
above. 

Although these above-mentioned papers constrained a 
wide variety of cosmological parameters, the compression 
was in all cases optimized only for one parameter, the 
overall power spectrum normalization. Although this pro- 
cedure can be improved as was described in Section p.6[ 
this does not appear to be causing substantial loss of ef- 
ficiency in the COBE case. Support for this conclusion is 
provided by Tegmark & Bunn (1995), where it is found that 
KL-compression optimized for measuring the normalization 
gives error bars on the spectral index n that are less than 
10% away from the minimal value that is obtained without 
any compression at all. It should also be noted that if one 
optimized the KL-compression for a different parameter 9, 
say 9 — n, then C would no longer be of the simple form 
of equation (^). Thus the signal-to-noise eigenmode treat- 
ment no longer applies, and the more general equation ( |29[ ) 
must be used instead. 

As has been pointed out by Gorski (1994), the rapid 
fall-off of the COBE window function We implies that vir- 
tually all the cosmological information in the COBE maps is 
contained in the multipoles I < 30. When implementing the 
KL-compression in practice, it is useful to take advantage 

t In fact, both Bond and Bunn independently reinvented the 
entire KL-method. 



of this fact by replacing the data set by the 957 spherical 
harmonic coefficients ae m for 2 < I < 30, since this reduces 
the size of the matrix to be diagonalized by about a factor 
4. This is an example of what we will call pre- compression, 
and we will return to this topic below, in Section |H| 

4-2.2 Signal-to-noise analysis of galaxy surveys 

The application of the signal-to-noise eigenmode method to 
galaxy surveys was pioneered by Vogeley (1995) and has 
since been further elaborated by Vogeley & Szalay (1996). 
Although these are both method papers, deferring actual 
parameter estimation to future work, the former explicitly 
evaluates the eigenmodes for the CfA survey and shows that 
there are about 10 3 modes with a signal-to-noise ratio ex- 
ceeding unity. 

In galaxy survey applications, the noise matrix N con- 
tains the contribution from Poisson shot noise due to the 
finite number of galaxies per unit volume, rather than on de- 
tector noise as in the CMB case. With galaxies, the question 
of what to use as the initial data vector x is not as simple as 
in the CMB case, since there is no obviously superior way 
to "pixelize" the observed density field. Vogeley & Szalay 
(1996) divide space into a large number of disjoint volumes 
( "cells" ) , and derive the resulting signal matrix S in the ap- 
proximation of ignoring redshift distortions. One alternative 
is using "fuzzy" cells as discussed by Tegmark (1995), for in- 
stance averages of the galaxy distribution where the weight 
functions are Gaussians centered at some grid of points. 

Another alternative, which has the advantage of mak- 
ing it easier to include the effect of redshift distortions, is 
to use the coefficients from an expansion in spherical har- 
monics and spherical Bessel functions as was done by HT95, 



Ballinger et al. (1995; hereafter BHT95) and in Section 4.1 
above. In BHT95 the real space power spectrum of density 
perturbations was left as a free function, to be evaluated 
in a stepwise fashion, along with the distortion parameter. 
In the case of only estimating the power, by fixing (5 to 
some constant value, the problem reduces to a signal-to- 
noise eigenmode problem, but with m free parameters — 
the power at m specified wavenumbers. Due to the mixing 
of wavenumbers by the <P and V matrices, these are gener- 
ally not independent and so we must apply the SVD method 
to construct orthogonal eigenmodes. This application will be 
discussed elsewhere. 

We conclude this section with a few additional com- 
ments on the above-mentioned "pixelization" issue. Com- 
pared to a CMB analysis, a galaxy survey analysis involves 
one more step, which is reducing the point data to the "ini- 
tial" data vector x. What is the relation between the first 
step (weighting galaxies) and the second step (weighting 
modes)? In Feldman, Kaiser & Peacock (1994) and Heavens 
& Taylor (1995), schemes for optimal weighting of galaxy 
data were presented. This weighting of the data (step 1) 
is prior to and complementary to the mode-weighting dis- 
cussed in this paper (step 2) . The data- weighting is designed 
to maximize the signal-to-noise of individual modes, and one 
might intuitively expect thus this to be a good mechanism 
for obtaining large eigenvalues A*,. The techniques of the 
present paper can then subsequently be used to trim the 
data in an optimal way. The data-weighting scheme alone 
does not guarantee that we get the best answer possible. 
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Figure 4. The three heavy lines show the smallest error bars 
attainable for the three parameters Q, n and r as a function of 
the number of modes used, i.e., the error bars obtained when using 
the separately optimized KL-modes for each parameter. The three 
thin lines show the error bars obtained when using the 500 SVD 
modes, illustrating that these contain essentially all the relevant 
information about all three parameters. 



This is because it maximizes the diagonal elements in the co- 
variance matrix, ignoring correlations between modes. This 
does therefore not guarantee that the complete Fisher infor- 
mation matrix will be optimal. In general, it does not seem 
tractable to devise a data-weighting scheme which optimizes 
the Fisher matrix directly. 

4.3 A multi-parameter CMB example 

Here we apply a KL-analysis to the problem of estimating 
the quadrupole normalization Q, the spectral index n and 
the reionization optical depth r from the 4-year COBE data 
(Bennett et al. 1996). As the initial data vector x, we use 
the temperatures in the N = 4016 pixels that l ie mo re than 





20° from the galactic plane. Just as in Section 4.2.1 



and C = S + N is given by equation (|?]). Computation 
of C,i thus reduces to computing the derivatives of the an- 
gular power spectrum Ce that occur in the sum, and we 
do this in practice using the software described in Seljak 
& Zaldarriaga (1996). The fiducial model has the standard 
CDM power spectrum shown in Figure |l| The three heavy 
lines in Figure ^ show the results of performing a separate 
KL-analysis for each parameter. In good agreement with 
the findings of Bunn & Sugiyama (1995), we see that vir- 
tually all information about Q is contained in the first 400 
modes. Similarly, we see that all essentially all the infor- 
mation about n and r is contained in the first 400 modes 
optimized for these parameters. Since the KL-modes for a 
parameter by construction give the smallest error bars, the 
curves corresponding to any other choice of modes would lie 
above these solid curves. For instance, if we used the KL- 
modes optimized for r to measure Q, the resulting curve 
would lie above the bottom one if we were to plot it in Fig- 
ure 0. Figure [l] provides a simple interpretation of this: since 
dCifdr « for I <C 10, the r-modes do not contain infor- 
mation about the lowest multipoles, since this information 



Figure 5. The error bars on the power spectrum normaliza- 
tion are shown for hypothetical COBE experiments with differ- 
ent noise levels. From top to bottom, they correspond to a noise 
enhancement factor 10, the real 4 year data, and noise reduction 
factors of 10, 100 and 1000. 



is useless for measuring r (even though it is important for 
measuring Q). 

We implement the SVD technique described in Sec- 



tion 3.£ by taking the 500 best modes for each parameter 



and SVD-compressing these down to 500 modes. The thin 
lines in Figure ^ show the resulting error bars. We see that 
although they lie above the minimal curves for k <C 500, 
they all "catch up" at k = 500. In other words, these 500 
modes retain virtually all the relevant information about Q, 
n and r, since using all of them gives virtually identical er- 
ror bars to those obtained when using the full N = 4016 
uncompressed data set. 



4.4 A low noise CMB example 

How effective will KL-compression be for future CMB data 
sets? Might it be the case that when noise levels are much 
lower, then all N KL-modes will have S/N3> 1, so that no 
compression is possible? Figure ^ shows that the answer to 
the second question is no. For this example, we have re- 
peated the above-mentioned COBE analysis for a range of 
different noise levels. With ten times less noise, the com- 
pression factor is seen to drop to about 4016/700 ~ 7 and 
another order of magnitude of noise reduction lowers the 
compression factor to about 4. The upcoming MAP experi- 
ment forecasts a noise level w -1 = 4-k<7 2 /N w 2.5 x 10 -15 , as 
compared to w^ 1 « 10" 12 for COBE and io _1 ?s5x 10~ 18 
for the planned COBRAS/SAMBA experiment. These are 
quadratic quantities, so taking square roots, we see that 
MAP and COBRAS/SAMBA reduce the COBE noise by 
factors around 20 and 450, respectively. Figure |E] shows that 
even with 1000 times less noise, a compression factor of 2 
is readily attainable. The explanation of this is oversam- 
pling. To avoid aliasing problems, the mean pixel separation 
must be smaller than the beam width by a substantial factor 
(the Shannon oversampling factor ~ 2.5). This redundancy 
remains even when the data set itself has excellent signal-to- 
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noise, so the KL-compression can take advantage of the fact 
that the pixel basis of the map is a sense overcomplete. Al- 
though future CMB experiments will of course have a much 
higher angular resolution than in our example, the oversam- 
pling factor will remain similar, so our conclusion about the 
usefulness of KL-compression will remain the same. 



5 ON GIANT DATA SETS: THE NEED FOR 
PRE-COMPRESSION 

Above we have shown that the KL-compression technique 
is in many cases very useful. In this section, we will discuss 
some of its limitations, as well as outline nonlinear exten- 
sions that may make the compression technique feasible even 
for the next generation of CMB missions and galaxy surveys. 

The number of pixels in a sky map from a next- 
generation CMB mission is likely to be around 10 7 . The 
Sloan Digital Sky Survey (SDSS) will measure the redshifts 
of 10 6 galaxies. Is it really feasible to apply KL-methods 
and likelihood analysis to data sets of this gargantuan size? 
We will argue that although an orthodox KL-analysis is not 
feasible, it appears as though the answer to the question 
may nonetheless be an encouraging yes if an extra nonlin- 
ear compression step is added to make the analysis doable 
in practice. Let us first note that despite the statistical or- 
thonormality property of a KL-compressed data set, the ma- 
trices that need to be inverted to find the ML-estimate are 
in general not diagonal. BCB 1 is only diagonal at the one 
point in parameter space where & = &o, i.e., the point 
corresponding to the parameter values that we assumed a 
priori. The ML-point generally lies elsewhere, and we need to 
find it by a numerical search in parameter space, so BCB* 
will generally not even be sparse. This forces us to resort to 
Cholesky decomposition. For the n ~ 4000 case of Tegmark 

6 Bunn (1995), this took about 10 minutes per matrix on 
a workstation. Since the time scales as n 3 and the storage 
as n 2 , this would require about 30 years of CPU time for 
n — 10 6 , and about one terabyte of RAM. Even allowing 
for the exponential rate at which computer performance is 
improving, such a brute force likelihood analysis does not 
appear feasible for a megapixel CMB map within the next 
ten years. 

The crucial question thus becomes how large a compres- 
sion factor we can expect to achieve. We will argue that a 
factor of 10 is just about all one can do with linear com- 
pression, but that quadratic "pre-compression" may be able 
to gain another factor of 1000 for a next generation CMB 
experiment. 



5.1 The limits of linear compression 

Consider the following simple one-parameter example: we 
are given n numbers Xi drawn from a Gaussian distribution 
with zero mean and an unknown variance 9 that we wish to 
estimate. Thus 



A9 = \ — 6 

n' 



(59) 



C 



0, 
91, 



(58) 



so when we solve equation ( |29[ ) , we find that all the eigenval- 
ues are identical: — 1/9. Thus if we compress by keeping 
only the first n' modes, the resulting error bar will be 



The fact that all eigenvalues are identical means that we 
would be quite stupid to throw away any modes at all, since 
they all contain the same amount of information. Simple as 
it is, this example nonetheless illustrates a difficulty with 
analyzing future CMB maps. Even in the best of all worlds, 
where we had complete sky coverage, we would encounter a 
problem of this type. To estimate a multipole Ci, we would 
be faced with (2£+l) coefficients ai m with zero mean and un- 
known variance, just as in the example above, which means 
that the KL-method would be of no use at all in compressing 
this data. 

In Bunn (1995) and in one of the above examples, it was 
found that the COBE data set can be compressed down to 
about 400 numbers without losing much information. Where 
does this magic number 400 come from? In the absence of 
a galaxy cut, the number would probably be around 600, 
since beam smearing has eliminated virtually all informa- 
tion about multipoles £ > 25, and there are about 600 ae m - 
coefficients with I < 25. The fact that the galaxy cut re- 
moves about one third of the sky then reduces the cosmolog- 
ical information by about a third. There is also a more direct 
way to see where the com pression factor 10 came from. As 
mentioned in Section 4.4 , pixelized maps are generally over- 



sampled by a factor 2.5 or more to avoid aliasing problems, 
which means that the mean pixel separation is considerable 
smaller than the beam width. Since the sky map is two- 
dimensional, the number of pixels per beam area is roughly 
the square of this number, of the order of 10. Future CMB 
missions will probably use similar oversampling rates, which 
means that a compression factor of 10 is probably the most 
we can hope for with linear compression only. 



5.2 Quadratic compression 

Fortunately, we are not limited to linear data compression. 
Let us compress the data set of our previous example into a 
single number defined by 



n 



(60) 



Let us assume that n is very large, say n > 100. Then y 
will be very close to Gaussian by the central limit theorem. 
This means that the mean is (y) = n9 and the variance is 
(y 2 ) — (y) 2 — 2n9 2 . Substituting this into Equation ( |l5| ) now 



A9 = 



n + 4 



(61) 



i.e., the theoretically minimal error bar of equation ( |59[ ) 
with n' = n. (The extra 4 in equation (|6l]), which might 
seem to indicate that one can attain smaller error bars than 
the theoretical minimum, should of course be ignored - 
it merely reflects the fact that the Gaussian approximation 
breaks down for small n.) 

In summary, we have found that whereas linear com- 
pression was powerless against our simple toy model, 
quadratic compression made mincemeat of the problem, con- 
densing all the information into a single number. This result 
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is hardly surprising, considering that the compressed data 
set y that we have defined in equation is in fact the 
ML-estimate of the variance. It nonetheless has far-reaching 
implications for the issue of how to compress future CMB 
maps. If we have complete sky coverage, and define the com- 
pressed data set as 
^ i 

ye = 21 + i l a «" l l 2 ' ( 62 ) 

m=-l 

then it is easy to show that the new Fisher information ma- 
trix will be identical to that of equation ([ll]), which involved 
using the entire data set. For an experiment with a FWHM 
beamwidth of 4 arcminutes, there is virtually no informa- 
tion on multipoles above f m ax = 3000, so this means that 
in the absence of a galaxy cut, we could compress the entire 
data set into 3000 numbers without losing any cosmological 
information at all. Roughly speaking, this works because 
the compression throws away all the phase information and 
keeps all spherically averaged amplitude information, and 
it is the latter that contains all the information about the 
power spectrum. (For non-Gaussian models such as ones in- 
volving topological defects, the power spectrum alone does 
of course not contain all the relevant information.) 

5.3 Real-world CMB maps: an outline of a 
compression recipe 

The discussion above already included the effects of pixel 
noise and beam smearing. Barring systematic errors, there 
are two additional complications that will inevitably arise 
when analyzing future high-quality CMB maps: 

• Foregrounds 

• Incomplete sky coverage 

Any attempts at foreground subtraction should of course be 
made before throwing away phase information, so that avail- 
able spatial templates of galactic dust, radio point sources 
etc. can be utilized. For a recent discussion of such sub- 
traction schemes, see e.g. Tegmark & Efstathiou (1996). 
The final result of the foreground treatment will almost cer- 
tainly be a map where some regions, notably near the galac- 
tic plane, have been discarded altogether. The resulting in- 
complete sky coverage degrades information in two different 
ways: 

• The sample variance increases. 

• The spectral resolution decreases. 

The former effect causes the variance in individual multipole 
estimates to grow as the inverse of the remaining sky area 
(Scott et al. 1994), and this increase in variance is automat- 
ically reflected in the final Fisher information matrix. The 
latter effect means that the ye defined by equation ( |62[ ) is no 
longer a good estimate of the multipole Ce . Rather, it is easy 
to show that (yt) will be a weighted average of all the mul- 
tipoles Ce- These weights, known as the window function, 
generally form quite a broad distribution around £, which 
means that the compressed data yt are effectively probing a 
smeared out version of the power spectrum. For a 20° galac- 
tic cut, this smearing is found to be around A£/£ ~ 25%, 
which clearly destroys information about features such as 
the exact location of the Doppler peaks. 



5.3.1 How to deal with incomplete sky coverage 

Fortunately, the problem of incomplete sky coverage can be 
for all practical purposes eliminated. As described in detail 
by Tegmark (1996, hereafter T96), it is possible to obtain 
much narrower window functions by simply replacing the 
spherical harmonic coefficients in equation ( |62| ) by expan- 
sion coefficients using a different set of functions, obtained 
by solving a certain eigenvalue problem. If is found that 
for a 20° galactic cut, the window functions widths can be 
brought down to A£ ~ 1 for all I, corresponding to a relative 
smearing of less than a percent at £ ~ 200, the scale of the 
first CDM Doppler peak. In other words, as long as none of 
the models between which we are trying to discriminate have 
any sharp spectral features causing the power spectrum to 
jump discontinuously from one multipole to the next, then 
virtually no information at all is lost in our quadratic com- 
pression. 

In general, smoothing only destroys information if there 
are features on the smoothing scale or below it. If the true 
power spectrum is similar to a CDM spectrum, it will typ- 
ically vary on scales A£ ~ Ce/ (dCe jdt) ~ 100, at least for 
angular scales £ > 50. In other words, even smoothing it with 
window functions with A£ as large as 10 would hardly de- 
stroy any information at all about this part of the power 
spectrum. The quadratic compression of T96 produces a 
spectral resolution of A£ ~ 1/8, where 9 is the smallest 
dimension of the patch of sky analyzed, in radians. In other 
words, we can allow ourselves even more leeway than the 
galaxy cut dictates without losing information. This can be 
used to save CPU-time in practice. Since we can achieve the 
spectral resolution A£ ~ 10 with 6° x 6° patches of sky, we 
can form compressed data sets y separately for a mosaic of 
such patches covering all the available sky, thus radically re- 
ducing the size of the matrices that we need to diagonalize 
for the T96-method, and then simply average these different 
estimates of the power spectrum. 

5.3.2 The final squeeze: KL-method 

Although the above-mentioned prescription will reduce the 
size of the data set dramatically, from perhaps 10 7 numbers 
to about 3000, there will still be considerable amounts of 
redundancy, since power spectrum estimates ye for neigh- 
boring ^-values will be correlated. Because of this, it is 
worthwhile to subject the new data set y to a regular KL- 
compression. We thus term the above-mentioned quadratic 
step pre- compression: it does by no means need to be opti- 
mal, and is simply done to reduce the data set down to a 
small enough size that it can be fed into the KL-machinery 
without practical difficulties. 

The values ye are Gaussian to an excellent approxima- 
tion for £ > 50, by the central limit theorem, since they are 
the sum of 2£ + 1 > 100 independent numbers. The remain- 
der, however, are not. Since the Gaussian property makes 
such a dramatic simplification both in the compression step 
and in the likelihood analysis itself, we therefore recommend 
discarding the y^-values with £ < 50 and replacing them by 
the 2500 spherical harmonic coefficients ae m with £ < 50 be- 
fore proceeding to the KL-compression step. This way, the 
entire data set y will be Gaussian. In addition, no infor- 
mation whatsoever has been lost about the largest angular 
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scales, where some models in fact do predict rather sharp 
spectral features which could render the quadratic compres- 
sion step inappropriate. 

In summary, we argue that the following prescription 
for analyzing future 10 7 pixel CMB map will be fairly close 
to optimal: 

(i) Foregrounds are subtracted making maximal use of 
multifrequency data and spatial templates obtained from 
other experiments. 

(ii) The most contaminated regions, for instance the 
galactic plane, are discarded altogether. 

(iii) The remaining sky is expanded in spherical harmon- 
ics up to I — 50 and the coefficients saved. 

(iv) This remaining sky is covered by a mosaic of over- 
lapping regions, whose diameter are of order 5° — 10°. 

(v) The angular power spectrum is estimated separately 
from each of these patches with the method of T96, up to 
I = 3000. 

(vi) All these power spectrum estimates are averaged. 

(vii) The first 50 multipole estimates are discarded, and 
replaced by the 2500 numbers from step (iii). 

(viii) The resulting data set (about 5500 numbers) is com- 
pressed with the KL-method. 

(ix) All cosmological parameters of interest are estimated 
jointly from this compressed data set with a likelihood anal- 
ysis. 

5.3.3 Open problems 

The above prescription for quadratic compression was nec- 
essarily rather sketchy, since a detailed treatment of these 
issues would go well beyond the scope of this paper. In- 
deed, the discussion above left several important questions 
unanswered, which are clearly worthwhile topics for further 
research. Here are two such examples: 

• How can information loss during pre-compression be 
minimized? Above we merely showed pre-compression to be 
lossless in the absence of noise (or with identical noise lev- 
els in all pixels). In the presence of noise, one would expect 
lossless compression to involve some form of inverse- variance 
pixel weighting, i.e., giving less weight to noisier pixels. On 
the other hand, pushing such noise weighting too far could 
broaden the window functions to the point where low spec- 
tral resolution led to irreversible information loss. 

• How is quadratic precompression best implemented nu- 
merically? For instance, are there particularly choices of 
shapes and sizes of the above-mentioned patches that sub- 
stantially facilitate the calculation of the final mean vector 
and correlation matrix? Is a direct analytic calculation of 
these quantities feasible (C involves terms quadratic in the 
power spectrum), or is it faster to resort to Monte Carlo 
techniques for this step? 

5.4 Real- world galaxy surveys 

Also for large future galaxy surveys, some form of precom- 
pression appears to be necessary before a KL-compression 
can be done. If we wish to retain all the information on clus- 
tering down to say 3ft" 1 Mpc scales in a 3D volume of typical 
dimension 300fe _1 Mpc, we clearly need to "pixelize" it into 



at least (300/3) 3 = 10 6 numbers. Since even the power spec- 
trum in the deeply non-linear regime contains valuable cos- 
mological information, it would not seem justified to simply 
ignore these small scales. 

Unfortunately, the galaxy survey problem is consider- 
ably more difficult than the corresponding CMB problem in 
that information about, for instance, redshift space distor- 
tions lies hidden not merely in the overall power spectrum, 
but also in the phases, in the differences between radial and 
transverse clustering patterns. 

As we have noted, transforming to weighted harmonic 
amplitudes is one way of reducing the number of data 
"points" without losing cosmological information For in- 
stance, in the case of the analysis of the 1.2Jy survey by 
HT95 and BHT95, only the first 1208 modes where used 
from the 2000 galaxies, with the limit on modes used being 
set by the survey size and the smallest scale still in the lin- 
ear regime. However, if the full range of available modes is 
desired, we might need higher order compression options. 

One way to implement quadratic compression, while re- 
taining the important phase information, is to transform to 
some coordinate basis which is orthogonal in redshift -space. 
We can then apply the quadratic compression to estimate 
a "power spectrum" in the transformed space, having lost 
none of the underlying phase information. Some progress 
along the lines of finding an orthonormal basis function 
for redshift space has been made by Hamilton & Culhane 
(1995). However, these methods still fail to adequately deal 
with shot-noise and the phase mixing introduced by a finite 
angular mask function. 

Thus the issue of whether one can do even better with 
galaxy surveys, while staying within the realms of numerical 
feasibility, remains a challenging open question. 



6 CONCLUSIONS 

We have given a comprehensive discussion of how to best 
estimate a set of cosmological parameters from a large data 
set, and concluded the following: 

• Well-known information-theoretical results roughly 
speaking state that a brute-force likelihood analysis using 
the entire data set gives the most accurate parameter deter- 
mination possible. 

• For computational reasons, this will be unfeasible for 
the next generation of CMB maps and galaxy surveys. 

• The solution is to use a good data compression scheme. 

• The optimal linear compression method is the 
Karhunen-Loeve method, of which the so-called signal-to- 
noise eigenmode method is a special case. 

• Although the standard KL-method applies only when 
estimating a single parameter, it can be generalized to the 
multi-parameter case by simply adding a step consisting of 
a singular value decomposition (SVD). 

• This SVD step also provides a simple way out of the 
Catch-22 situation that one needs to specify the parameter 
values before one has measured them. 

• The KL-method produces a compression factor ~ 10 
for typically sampled CMB maps, and also for the redshift 
space distortion analysis of Heavens & Taylor (1995). 

• However, this is not enough to handle a high-resolution 
next generation CMB map. 
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• Fortunately, it appears as though this can be remedied 
by adding a quadratic pre-compression step without sub- 
stantial information loss. 

Cosmology, which used to be a data-starved science, is now 
experiencing a formidable explosion of data in the form of 
both CMB maps and galaxy redshift surveys. Around the 
turn of the millennium, we will probably be equipped with 
data sets so rich in information that most of the classical 
cosmological parameters can — in principle — be deter- 
mined with accuracies of a few percent or better. Whether 
this accuracy will be attainable also in practice depends cru- 
cially on what data-analysis methods are available. We have 
argued that the prospects of achieving this accuracy goal 
are quite promising, especially on the CMB side (which is 
slightly simpler), by using a multi-parameter generalization 
of the Karhunen-Loeve method together with a quadratic 
pre-compression scheme. However, much work remains to be 



done 
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